/*
  Copyright (C) 2011, Kyungjoo Kim
  
  This file is part of LINAL (LINear ALgebra)
  
  All rights reserved.

  This library is free software; you can redistribute it and/or
  modify it under the terms of the GNU Lesser General Public
  License as published by the Free Software Foundation; either
  version 2.1 of the License, or (at your option) any later version.
  
  This library is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  Lesser General Public License for more details.
  
  You should have received a copy of the GNU Lesser General Public
  License along with this library; if not, write to the Free Software
  Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA
  Also add information on how to contact you by electronic and paper mail.
  
  
  Kyungjoo Kim
  iamkyungjoo@gmail.com
*/
#include "linal/common.hxx"
#include "linal/util.hxx"

#include "linal/format.hxx"

#include "linal/operation.hxx"

#include "linal/device/host/operation.hxx"

namespace linal {
  // FLA_QR2_UT_inc_blk_var1
  int qr2(Hier_ U, Hier_ D, Hier_ T) {
    FLA_Obj UTL,   UTR,      U00, U01, U02, 
            UBL,   UBR,      U10, U11, U12,
                             U20, U21, U22;

    FLA_Obj DL,    DR,       D0,  D1,  D2;
    FLA_Obj TL,    TR,       T0,  T1,  W12;

    FLA_Part_2x2( ~U,    &UTL, &UTR,
                         &UBL, &UBR,     0, 0, FLA_TL );
    FLA_Part_1x2( ~D,    &DL,  &DR,      0, FLA_LEFT );
    FLA_Part_1x2( ~T,    &TL,  &TR,      0, FLA_LEFT );

    while ( FLA_Obj_min_dim( UBR ) > 0 ){
      FLA_Repart_2x2_to_3x3( UTL,  UTR,       &U00,  &U01, &U02,
                                              &U10,  &U11, &U12,
                             UBL,  UBR,       &U20,  &U21, &U22,
                             1, 1, FLA_BR );
      FLA_Repart_1x2_to_1x3( DL,   DR,        &D0,  &D1, &D2,
                             1, FLA_RIGHT );
      FLA_Repart_1x2_to_1x3( TL,   TR,        &T0,  &T1, &W12,
                             1, FLA_RIGHT );
      //------------------------------------------------------------
      Hier_ U_11(U11), D_1(D1), T_1(T1);
      FLA_QR2_UT_internal( U11,
                           D1, T1, 
                           fla_qr2ut_cntl_leaf );

      if ( FLA_Obj_width( U12 ) > 0 ) {
        Hier_ U_12(U12), W_12(W12), D_2(D2);
        FLA_Copy( ~U_12.flat(), ~W_12.flat() );
        FLA_Gemm( FLA_CONJ_TRANSPOSE, FLA_NO_TRANSPOSE,
                  FLA_ONE, ((FLA_Obj*)D_1.get_buffer())[0],
                  ~D_2.flat(), FLA_ONE, ~W_12.flat() );

        FLA_Trsm( FLA_LEFT, FLA_UPPER_TRIANGULAR, 
                  FLA_CONJ_TRANSPOSE, FLA_NONUNIT_DIAG, 
                  FLA_ONE, ((FLA_Obj*)T_1.get_buffer())[0],
                  ~W_12.flat() );
          
        FLA_Axpy( FLA_MINUS_ONE, ~W_12.flat(), ~U_12.flat() );

        FLA_Gemm( FLA_NO_TRANSPOSE, FLA_NO_TRANSPOSE, 
                  FLA_MINUS_ONE, ((FLA_Obj*)D_1.get_buffer())[0],
                  ~W_12.flat(), FLA_ONE, ~D_2.flat() );
      }
      //------------------------------------------------------------
      FLA_Cont_with_3x3_to_2x2( &UTL,  &UTR,       U00, U01,  U02,
                                                   U10, U11,  U12,
                                &UBL,  &UBR,       U20, U21,  U22,
                                FLA_TL );
      FLA_Cont_with_1x3_to_1x2( &DL,   &DR,        D0, D1,  D2,
                                FLA_LEFT );
      FLA_Cont_with_1x3_to_1x2( &TL,   &TR,        T0, T1,  W12,
                                FLA_LEFT );
    }

    return true;
  }
}
